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Abstract 

We investigate the extension of the Prokof 'ev-Svistunov worm algorithm 
to Wilson lattice fermions in an external scalar field. We effectively simulate 
by Monte Carlo the graphs contributing to the hopping expansion of the two- 
point function on a finite lattice to arbitrary order. Tests are conducted for 
a constant background field i. e. free fermions at some mass. For the method 
introduced here this is expected to be a representative case. Its advantage 
is that we know the exact answers and can thus make stringent tests on 
the numerics. The approach is formulated in both two and three space-time 
dimensions. In D = 2 Wilson fermions enjoy special positivity properties 
and the simulation is similarly efficient as in the Ising model. In D = 3 
the method also works at sufficiently large mass, but there is a hard sign 
problem in the present formulation hindering us to take the continuum limit. 
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1 Introduction 



To give a definition of quantum field theories beyond perturbation theory (Feyn- 
man diagrams) they normally have to be regularized by replacing space and usually 
also time by a lattice. Then the functional integral becomes a well-defined object 
and thus amenable to numerical methods, usually in the form of stochastically 
sampling lattice field configurations by Monte Carlo methods. Most theories of 
interest contain fermion fields which lead to some or all integration variables being 
anticommuting Grassmann 'numbers'. In the standard approach these integrations 
- possibly after introducing additional Bose fields — are Gaussian and are per- 
formed exactly. The result is an effective action of the bosonic fields alone which 
become coupled non-locally. The known Monte Carlo techniques to practically 
simulate such systems are mostly based on molecular dynamics and the hybrid 
Monte Carlo idea (HMC) [Tj, [2J. These methods have been improved and opti- 
mized rather successfully over the years by a very large effort of many members of 
the lattice community. On the other hand practitioners know that once fermions 
are decoupled then HMC for locally coupled Bose fields is not a very efficient algo- 
rithm compared to alternatives like over-relaxation which are then available, not to 
mention the (unfortunately few) cases where cluster methods can be applied. This 
implies a large penalty for fermions even in cases where their effects are only small. 
At small fermion masses the fermionic forces in HMC tend to grow and the step 
size of the molecular dynamics trajectories has to be taken small enough. With 
this quasi-continuous evolution one then has to be cautious about possible long 
autocorrelations. After all, to the best of the author's knowledge, the ergodicity 
of HMC has not been formally proven. 

Maybe for the aforementioned reasons among others some part of the commu- 
nity has remained motivated to look for radically different approaches. A rather 
natural idea is to look for a representation of fermions as some sort of 'sum over 
configurations' more similar to the bosons. One of the pioneering papers develop- 
ing such ideas is [3J. There as in numerous succeeding attempts one starts from an 
operator formulation of fermions and inserts intermediate states in the occupation 
number basis between factors of the transfer matrix. In this way occupied sites 
map out an ensemble of 'world-lines' or a gas of loops of fermions on the lattice. 
Often the amplitudes that arise oscillate in sign with the danger of leading to an 
unmanageable signal to noise ratio, the infamous fermionic sign problem. The 
inclusion of gauge fields in this approach poses additional problems. 

A somewhat different approach was successful for — but also restricted to - 
QCD at infinite gauge coupling, (3 — 0, [3]. In the Euclidean path integral with 
staggered fermions but no gauge plaquette term the group valued gauge fields can 
be integrated out first. The resulting model of locally paired even Grassmann 
elements has contributions that can be viewed as a statistical system of explicitly 
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color-neutral mesons (dimers) and baryon loops. Both these systems and the 
world-line gas a la [3] are difficult to simulate efficiently by local methods due to 
constraints which conflict with local deformations of the configurations. In some 
cases efficient nonlocal updates could be devised [5], [6]. 

A 'more Euclidean' version of the idea was proposed in [7j. These authors 
started from the determinant of the integrated-out staggered fermions and tried 
to stochastically generate its expansion into cycles. However the restriction to 
local updates and the sign problem, even for D = 2 in this case, have limited the 
use of the method. 

In [S] (see also [§]) a loop or world-line representation was proposed for the 
partition function of standard two-dimensional Euclidean Wilson fermions in an 
external scalar field. They were mapped on a certain 8-vertex model. Based 
on it the Gross-Neveu model was simulated with local updates in [10J. In [TT] 
the same representation was re-derived directly from the Grassmann integral for 
charge self-conjugate (Majorana) Wilson fermions. The mapping between Wilson 
fermions and a loop-gas could in addition be made precise also for a finite torus 
with (anti)periodic boundary conditions. A cluster algorithm for the loop-gas was 
developed in [TT] which produces almost uncorrelated loop configurations at low 
cost. In the sequel Willi Rath and the author have tried to compute correlations 
based on these configurations [12] . The only solution we have found so- far proceeds 
via the numerical generation of the scalar cr-field that usually factorizes the Gross- 
Neveu interaction. Then the close to singular Dirac operator in this random scalar 
field has to be inverted and the CPU time ends up being spent in a very similar 
fashion as in HMC 

In this paper, as an alternative approach, we adopt the 'worm' algorithm of 
Prokof 'ev and Svistunov (PS) [13] to lattice fermions of the Wilson type. To this 
end we build on the study of the PS algorithm for the Ising model carried out 
as a preparation in [14J. While there the (untruncated) strong coupling expan- 
sion is sampled, the fermion loop-gas corresponds to the quite similar hopping 
expansion^- We here extend the loop-gas formulation of fermions on a torus in 
two ways. We generalize [TT] to including two spinor field insertions at arbitrary 
lattice sites. It turns out that the PS algorithm is ideally suited to keep track of 
the non-local amplitudes involved due to Fermi statistics. The second non-trivial 
extension takes this construction to Majorana fermions in three Euclidean dimen- 
sions. While we understand why the fermionic sign problem mentioned before is 
absent in two dimensions if the system size is large in correlation lengths, the full 
problem has to be confronted in three dimensions. We indeed find for free fermions 
that are implemented numerically in this study, that the PS algorithm for D = 2 
is similarly efficient as in the Ising model. While clearly correct in principle also in 

1 Because of this strong similarity, we put this paper into one series with |14) . 



3 



D = 3 it fails numerically with the present technique when the continuum limit is 
approached. We nonetheless find the three dimensional loop representation theo- 
retically quite interesting. We think that the free Majorana fermion in D = 3 is an 
excellent study ground for more clever techniques, for instance cluster improved 
observables, to still overcome the sign problem, perhaps along the lines of |15j . 

The organization of this paper is as follows. In the next section we set up our 
notation for the lattice fermions discussed followed by section 3 introducing dimers 
that label all possible hopping graphs needed for the PS simulation. Tools for the 
simulations are described in 4. In section 5. we define the kind of observables on 
the loop ensemble that allow to make contact with fermionic two-point functions 
followed in 6. by the description of numerical results. We end on 7. conclusions 
including a brief outline how interaction can be added. In two appendices we 
collect the free fermion results used as benchmarks and a geometrical discussion 
of the fermionic phase factors arising for each closed loop. 

2 Majorana- Wilson lattice fermions 

We start from a standard Wilson-Dirac fermion with the action 



We consider a D- dimensional standard hypercubic lattice with spacing a in all di- 
rections and either periodic or antiperiodic boundary conditions for each direction 
over the respective periodicity length L M . The boundary conditions are coded into 
a vector with components 0,1 by the condition 



and similarly for ip, and ft is a unit vector in the positive /z direction. 

Unless stated otherwise, the mass m is assumed to be a real x-dependent 
periodic external field m(x) here. By later integrating over it with a suitable 
weight one can, starting from this building block, arrive at interacting theories like 
the Gross- Neveu model. The operators d,d*,d are the usual forward, backward, 
and symmetrized nearest neighbor differences. The set {7^, \i — 0, 1, . . . , D — 1} 
are hermitean Euclidean Dirac matrices. From here on we shall restrict ourselves 
to the space-time dimensions D = 2, 3 with 2x2 7-matrices in both cases. The 
Wilson term suppresses the doublers and from here on we set its coefficient to the 
convenient value r = 1. 

The action ([T]) is invariant under charge conjugation for any m(x). It is hence 
both possible and natural to split the fermion into two neutral Majorana compo- 




(1) 



X 




(2) 
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nents by setting 



^ = -L(£ 1 + ^ 2 ), ^ = -^=(£- i £)C (3) 



with the charge conjugation matrix C obeying 

c^c- 1 = - 7 j = - 7 ;, c = -c T . (4) 

Inserting this into (pQ) we find two identical contributions for £ 1)2 . In our Majorana 
reduction we consider only one such component in the following 

S = \a D Y, eC(lA +m- \ad*d% (5) 

X 

Note that the matrix in this quadratic form is antisymmetric. By collecting diag- 
onal and neighbor terms we can rewrite this action as 

s = I E( D + ™)e T m(x) - J2 C(x)cp(mx + # (6) 

X X,fl 

where we now have adopted lattice units (a = 1) and have introduced projectors 

P(n)= 1 -(l-n flltl ) (n 2 = l) (7) 

for each lattice direction (n = ±/t). Note that the hopping term of a Majorana 
fermion is a function of the unoriented link because of the identity 

eix)cp(fi)^ + a) = + p)cp{-m{ x ). (8) 

For D = 2 the form coincides with the starting point of [11J. 
To continue we introduce the shorthand notation 

z=ec (9) 

We emphasize that for the Majorana fermion this depends on £ while were 
independent Grassmann integration variables. The partition function is given by 

= J D£e- S = Pf [C( 7 A + m - \d*d)} (10) 

where the Gaussian integral over Majorana fields has led to a PfaffiarE of the 
antisymmetric matrix. The result depends on the boundary conditions, of course, 
which is exhibited for but left implicit on the right hand side. 



2 The order of factors in D£ is assumed to be such that this is true without an (irrelevant) 
extra sign. 



5 



In a straight-forward generalization of [M] we now extend our study to include 

Z (£ \u,v) = J Die- S i{u)l{v), (11) 
which is a matrix in spin space. It is closely related to the two point functional 



GU.U-.n,)- (i.U)V!,)) = J % V) . (12) 



As we are considering bilinear fermions in an external field m(x) the propagator 
can also be obtained as the solution of a system of linear equations 

(7/A + 171 ~ ^ d *d)G(x, y; m) = 8 x>y x l spin (13) 

where the Dirac operator acts on x. 

For constant m such an evaluation can proceed by Fourier expansion and will 
serve us as a check below. Otherwise the Pfaffian is a problem similar to the fermion 
determinant and methods like HMC are suitable at least for an even number of 
flavors [16], [TTj . Our objective here is however to develop a simulation method 
alternative to this approach. 

We end this section with the remark that, in contrast to the Ising model, 
not equal to the ordinary partition function. Instead one may show 
that for any m 

r)7 {£) 

zW ^> = a^) xl - (14 > 

holds. To derive this relation we use that the space of antisymmetric 2x2 matrices 
is only one-dimensional, given by multiples of the second Pauli matrix. Hence in 
;) the integral containing £,(x)l;(x) T must be proportional to C 1 . 



3 Dimer form of Majorana fermions 

We here derive the loop-gas form of the fermion correlation function and partition 
function. In principle this may be achieved by using theorems for the expansion 
of the Pfaffian together with the sparseness of the matrix introduced before. This 
would parallel the approach in [7]. Instead we shall extensively manipulate the 
representation by a Grassmann integral. This is physically more transparent and 
may be seen as deriving the required expansion formulas 'on the fly' as they are 
needed. The power of Grassmann numbers for such purposes was emphasized 
before in [T8] . 

3 The dependence of G on the boundary conditions e is left implicit. 
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3.1 General structure 



We start from the factorized form 



J z 

with the short hand 



e 2 



¥>(*)€(*)€(*) 



l=(xy) 



i{x)P(f^)i(y) 



£{ U )Z(v) (15) 



(p(x) = D + m(x). (16) 

Because each P is a one-dimensional projector and due to the Grassmann nature 
of £ there are only two termg^l in the expansion of each link-factor. It can thus be 
'dimerized' 



J2 [C(x)P(^x)C(y)] kl 



(17) 



h=o,i 



leading to 



Z^{ 



u, v 



{h} J 



e 2 



JJ [£(x)P(^)£(y)]' 

l={xy) 

For each configuration {/;;;} we say that links with k\ = 1 carry an (active) dimer. 
Associated with each site we have only two Grassmann variables integrated over. 
This implies numerous constraints on contributing dimer configurations: 

• at sites x {u, v} there can only be either or 2 dimers adjacent 

• if u v, at these two sites there must be exactly 1 dimer touching 

• at u = v there can be no dimer touching. 

In the simpler case of Zq £) the analogous expansion requires or 2 dimers around 
all sites. Any dimer configuration that obeys these conditions and contributes to 
(u,v) or to Zq we call admissible. 

As a consequence of these constraints dimers in admissible configurations have 
to form chains. These can never backtrack, intersect or overlap. For u ^ v there 
must be exactly one chain or string connecting u and v that we call a. Apart from 
it all other chains must form a number of closed loops Xj. For contributions to 
Zq there is no string but only loops and the same is true for Z^(u,u) with the 
additional requirement that no loop passes through u. 

It is to be emphasized that the string (if present) and the loops including their 
number is a unique one-to-one representation of an admissible dimer configuration, 

{h} I admissablc < > & U {Xj, j = 1, 2, . . . , A^ A }, (19) 

where both the string and the set of loops can also be empty. 



'While many of the previous steps go through also for D = 4, we will need a third term here. 
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3.2 Amplitudes 



To evaluate a contribution for fixed admissible {h}, we can reorder freely all Grass- 
mann bilinears like hopping terms of links with ki = 1, the local terms and the 
integration measure with two spin components d^d^ for each site. In this way 
the whole Grassmann integral is factorized into one site integrals that are carried 
out with the formula 

/ = Ispin, (20) 

where we have chosen C 12 = 1 = — C 2 i. At each monomer site — a site with 
no dimer adjacent — the integrations are saturated by a site-factor in (|18|) and 
contribute a factor (p(z) = D + m(z). 

Next we consider a and define |er| to be the number of dimers one has to cross to 
walk from u to v. During the walk one encounters a sequence of sites Sj separated 
by lattice unit vectors rij, 

string a {u = s ,s 1 ,s 2 , ... ,s\ ff \ = v}, n i+1 = s i+1 - s { . (21) 

We notice that we may use © to our convenience along the path. Then, after 
carrying out the integrations belonging to all sites s i} there emerges a product of 
projectors 

V{a)=P{n 1 )P{n 2 )---P{n W \). (22) 

Each closed loop Xj can be labeled in exactly the same way except that now 
so = S|a| holds (suppressing for the moment the loop index j). If we denote the 
sequence of unit shifts now by {m 1 , m 2 , . . . , rri\ X \} , then closed loops contribute a 
scalar factor 

w(Xj) = - tr[P(m 1 )P(m 2 ) • • • P(m w )]. (23) 

The minus sign here is the usual one coming from closed fermion loops. Technically 
speaking, upon closing the trace, one pair of £, £ appears in the 'wrong' order. The 
cyclicity of the trace immediately implies, that w(Xj) is independent of where we 
start with mi along the loop. In addition one may use 

P{n) T = CPi-n)^ 1 (24) 

to show also independence of the direction chosen to traverse the loop. Hence 
w(Xj) is truly a function of the unoriented loop only. 

3.3 Evaluation of spin factors 

Using a bra-ket notation in spin space we write the Wilson projectors as 

P(m) = \m)(m\, (m\m) = 1, m = ±/2 = ±6, ±1[, ±2 if D = 3]. (25) 
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Now the loop factor is composed of scalar products 

w (^j) = -{mi\m 2 ) (m 2 |m 3 ) • • • (m|A|-i|ra|A|) (m\\\\mi) , (26) 

i.e. factors associated with the sites connecting pairs of links met along the loop. 
The modulus of the individual factors is one where successive mi coincide (straight 
sections) and l/v2 = cos(7r/4) between orthogonal links (corners). This is easily 
seen from a simple example 

|(6|i)| 2 = <6|P(i)|6> = (6| 7o P(i) 7 o|6) = (6|P(-i)|6) = \ (27) 

and similarly for any other orthogonal pair. With 7r/4 we see the typical half-angle 
appearing with spinors. We thus find 

w{\) = 2-^/ 2 <KA,), (28) 

where C(Xj) is the number of corners around the loop and 4>(\j) is a phase. 

We now discuss a rather direct way to compute <p{\j)- While it gives not 
much geometric insight into its meaning, this derivation will lend itself to a very 
direct algorithmic implementation. In appendix [B] an alternative more geometrical 
analysis is presented. 

We fix the ambiguous phases of | ± (x) in a definite way, knowing that w is 
independent of this convention. We start from |0) with an arbitrary phase. Next 
one may demand a maximal number of five real positive phases 

(±i|6) = <-6|i) = (±2|6) = -L (29) 

This exhausts the free choices and the remaining phases, one in D = 2 and addi- 
tional six in D = 3 can be evaluated. One possible way to do so is to construct all 
eigenvectors starting from |0) with the help of the projectors: 

|±i)= v / 2P(±i)|6), (30) 

|-6)=v / 2P(-6)|l) = -7i|6) (31) 

and 

| ±2)=v / 2P(±2)|6). (32) 

The implied phases can now be computed by just using the Dirac algebra and 
they are collected in table [TJ In two dimensions only the upper left 4x4 block is 
relevant, where all phases are real, they belong to Z(2). In this case the sign for 
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Table 1: Phases in all possible scalar products between eigenspinors. 



each loop acquires a simple geometrical interpretation which will be discussed in 
section [5l 

In evaluating the phases involving the third dimension we have assumed that 
7i72 = — ^7o holds and set 

1 + i i » . . 

z=^— = e l *. 33 

We cannot eliminate the complex phase factors by re-defining the phases of | ±2). 
There is another inequivalent irreducible Dirac representation in D = 3 with the 
complex conjugate phases (for example from 7 M — > —7^). This means that the 
parity reflected loop has the opposite phase which is hence a pseudoscalar. The 
phase factors for 3-dimensional Wilson fermions are in Z(8) and all values are 
actually assumed for relatively short loops if non-planar ones are included. The 
group Z(8) is related to the rotation group being reduced to the lattice symmetries, 
see appendix [Bl Simple examples for loops with complex phases are shown in figure 
HJ They were extracted from a Monte Carlo simulation (see below) be 'tagging' 
the phase of configurations and plotting them. 



Figure 1: Closed fermion loops in D = 3 with phase exp(i7r/4) (left) and exp(i7r/2) 
(right). 

The spin factor for the open string from u to v is given by 

V(a) = |ni)(ni|n 2 ) ■ ■ • (n H _i|n H ) (n\ \\. (34) 
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For the leftmost ket and the rightmost bra we introduce the notation 



|m) = \n(u)}, (n\ a \\ = (n(v)\ (u ^ v) (35) 

such that n(u) is the unit vector pointing out of u in the direction of the unique 
adjacent dimer k\ = 1 while n(v) is the corresponding unit vector pointing toward 
v. Note that in principle we should write n(u; k) and n(v ; k) and both vectors are 
undefined if u — v holds. The scalar factors have again a modulus for each 
corner and a phase <fi{a) that may be constructed from table 

V{&) = 2-°^<p(a) |ti(m))<7i(i;)|. (36) 

We re-emphasize that all objects discussed above including the number of cor- 
ners, the string and loop decomposition and the various phases are (nonlocal) 
functions of the k\ in an admissible configuration. It would however clutter our 
notation too much to always exhibit this explicitly. 

3.4 Boundary conditions 

More phase factors can arise from boundary conditions if loops or the string winds 
around the torus in antiperiodic directions an odd number of times. We adopt the 
convention to label the points on the torus by coordinates = 0, 1, . . . , L M — 1 
and distinguish a (D — l)-dimensional sheet of 'boundary' linkqfl for each direction 
as follows: 

/ is a boundary link in direction /i <-> I = (x,x + p) with x M = L M — 1. (37) 

For the string a and for each loop Xj we introduce parities e^{a) and e At (Aj) defined 
by 

, s f 1 if a contains an odd number of u — boundary links 
6 > } = \ else (38) 

and for the closed loops e At (A J -) is completely analogous. The overall sign from the 
boundary conditions is now given by 

sign = (-l)^ with e lt = e (t {a) + Y^^(X j ) (mod2) (39) 

3=1 

with the scalar product of the D- vectors e M and e M in the exponent. 

All contributions to the amplitude of an admissible dimer configuration have 
now been identified and will be combined in the next section to represent fermionic 
quantities in an ensemble summing and ultimately sampling such configurations. 



3 Of course, the torus has no boundary, hence the quotes. 



11 



4 Dimer partition function and worm algorithm 
4.1 Dimer partition function 

Next we formally write down the characteristic function Q(k; u, v) that is unity for 
admissible configurations and zero for all others. In principle it has been defined 
before in words. As a building block we use 

d(k; x)=J2 k i ( 40 ) 

l,dl3x 

which counts the number of dimers adjacent at x. Then we define 

Q(k;u,v) = 6 d {k;u),ibd{k;v),i ] f {Sd(k;x),o + 8d(k;x),2) for u ^ v (41) 

x£{u,v} 

Q(k;u,v) = Y[(5 d {k;x),o + S d {k;x),2) for u = v. (42) 

X 

Note that the constraint enforced for u = v here is the one for contributions to 
Zq rather than the more restrictive one for contributions to Z^(u,u). We now 
consider the following partition function 

Z= Y] ^hA 2 ~c/2 TT ( x). (43) 
^ p(u,v) 11 

u,v,{ki} ry ' x,d(k;x)=0 

Here C is the total number of corners 

C = C(a) + J £C(X j ), (44) 
i=i 

p is an arbitrary symmetric strictly positive lattice-periodic function similar as 
in [H]. The product is the weight from all monomer sites and we here restrict 
ourselves to 

tp(x) = D - m(x) > (45) 

guaranteeing the positivity of the overall weight. Expectation values of observables 
A(k;u,v) in this ensemble are defined by 

(W) = 7 E A(t; ¥ )fe!2^ J] (46) 
Z ^-^ P[u,v) ±J - 

u,v,{h} rv ' ' x,d(k;x)=0 

Observables related to the Majorana fermions discussed before can be written 
as ratios of such expectation values. This will be discussed in the next section 
after introducing the simulation algorithm for (|43|) . 
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4.2 Prokof 'ev-Svistunov worm algorithm 



The simulation of the dimer ensemble can be carried out with the worm algorithm 
of PS [13]. It is very similar to the algorithm described in [14J and we can be 
brief here about details. The main difference to the Ising case is that more than 2 
dimers cannot touch and that there is a weie ht l/y/2 for corners which induces a 
kind of stiffness (tendency to be straight) of the chains. 

We briefly pause here to comment on the notation of hopping parameter ex- 
pansion in the title of the paper. The factors <p(x) appearing for monomers could 
be rescaled to unity by absorbing them into £(x) early on. Then each dimer k/ xy \ 
would be accompanied by a factor [Lp( y x)f(y)]~ 1 ^ 2 . For constant m this would 
equal [D + m] _1 = 2k with the conventional hopping parameter k. Thus 2k is 
the strict analog of the Ising strong coupling expansion parameter tan/3 in |14j . 
We prefer however to stay with the unrescaled form which is advantageous for the 
introduction of interaction via m(x). 

An update microstep of the PS algorithm is now a succession of steps I and 
II applied to admissible configurations. In step I we make a Metropolis decision 
on a proposal where we pick one of the 2D nearest neighbors of v with equal 
probability and call it v' and the connecting link I. The proposed move changes 
v — > v' flipping at the same time fcj — > k[ — 1 — hi. It brings us from the global 
configuration k to configuration k' (differing at exactly one link). Note that k! 
may be not admissible, in which case the move will be rejected. We first form an 
auxiliary quantity q, the ratio of amplitudes after and before the move. We have 
to distinguish a number of cases and collect values of q in table [2j We have written 
n(u') as a shorthand for what should be n(u] k') after the move, although we did 
not move u here. The allowed moves are illustrated in figure [21 Parts a), b), c) 
refer to lines 1, 3, 5 of the table. The lines below those refer to the reverse changes. 
They correspond to the same graphs read from right to left with the arrow reversed 
and interchanged v <-» v'. The directions of active dimers not participating in the 
present update are examples only and can also point differently. 
In all other cases not covered here q is set to zero. In lines 3 and 4 of the table a 
sign is included for changing the number of fermion loops N\ by one. Finally the 
modulus of q is used in the acceptance probability 



while the phase changes will be considered in the next section. 

The type II move is as follows. If we encounter a configuration u, v, {k{\ with 
u = v we 'kick' u = v together to a randomly chosen other lattice site with 
unchanged {k{\ with the probability < p^ick ^ 1- For u ^ v we do nothing in 
this step, which is the dominant case. A difference with respect to the Ising case 




(47) 
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Figure 2: Pictorial representation of elementary moves in the PS algorithm. Solid 
lines are active dimers (value one). 

is that while there also pkick = (absence of step II) yields an ergodic algorithm 
|19j . this is not so here. For the fermions the jumps are required to move between 
different connected components. 

Moving only v together with steps II constitutes a correct Monte Carlo al- 
gorithm. We nevertheless found it advantageous to also move a in a completely 
analogous fashion. We thus now call the sequence l u — II — 1„ — II a microstep and 
call N x /2 mircosteps an iteration if we have N x lattice sites. 

We found the choice of pkick £ [0-3, 1] not critical and use pkick = 0.7 in the 
following after a few quick experiments. 



d(k; v) 
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Table 2: Entries in the first four columns specify the condition for possible moves, 
under which the amplitude gets multiplied by q (not including possible signs from 
antiperiodic boundary conditions). 
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5 Fermionic phase and spin factors 



5.1 Formulae for both D = 2 and D = 3 

The configurations of the dimer ensemble just discussed correspond to the set of 
graphs of the hopping parameter expansion of Majorana fermions. Each admissible 
configuration contributes to Zq or to Z^{u,v) with a certain amplitude and in 
the second case also with a spin matrix. The moduli of the amplitudes have been 
incorporated into the generation of configurations. The phases of the amplitudes 
and the spin matrices will be taken into account now as observables evaluated as 
in flU}. 

We first combine all phases discussed in section |3~B1 and I3~4l into the total phase 

*W(*) = (-1)»II^). (48) 

J=l 

For u = v we set 0(cr) = 1, there is no string, only loops. 

As mentioned before for D = 2 the phase $ is just a sign while for D = 3 it is 
an element of Z(8). We now have the connection 

VZ { £) = Z((p(u,u)5 u ^\k))) (49) 

with the volume 

D-l 

V = J] L M . (50) 

If we define a spin matrix 

, > / \n(u))(n(v)\ if u^v , , 

S{k;u,v) = < ' " \ y , j (51) 

the cases with insertions may be uniformly written as 

Z^\x, y) = p(x, y)Z((5 U) J V;y ^\k)S(k; u, v))) . (52) 

At coinciding arguments x = y the Grassmann integrations are saturated by the 
insertion alone, which requires a monomer site with its usual weight factor to be 
canceled. In this case one could in principle also relax the constraint in (|52p to 
obtain 

VZ {£ \x,x) = Z ((p(u,u)5 U:V $ {£) (k)S(k;x,x))) . (53) 
For the fermion correlation function the connection is 

_ p(x,y)((5 u J v ^(k)S(k;u,v))) 
G{X ^ m} - (1/V)({p(u,u)8 u , v mk))) ■ { } 
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An alternative derivation of the result for coinciding points starts from the obser- 
vation 



d 




1 {{p{u,u)5 u ^ £) (k)5 d ^ k - x)fi 



» 



(55) 



dm(x) 



and uses (H9l) and 

From the structure of the contributions in the dimer ensemble we may conclude 
that the right hand side of is rational in the external field m(x). The denomi- 
nator has total degree V, the numerator V — 1 — |er| m i n . Here |o~| m i n is the minimal 
number of links to connect x and y by a string. The degree in each individual 
m(x) is only linear both in the numerator and denominator. 

The above formulae simplify if translation invariance holds, m(x) = m, G(x, y; m) 
— > G(x — y), where we also restrict p(x, y) = p(x — y) and normalize p(0) = 1. We 
then find 



We recognize close similarities with the Ising correlation in [H] with the novelty 
of averaging the phasefactor and the spin matrices. Note that the delta function 
8& in the numerator needs to have the same antiperiodicity as the fields £. 
From ( H9i) one may now trivially read off that 



allows to measure the change in free energy for different boundary conditions. 

This type of quantity is theoretically nice, because it is expected to possess a 
continuum limit in a finite volume. For the massless all-periodic case = 0, the 
partition function vanishes at m = because the matrix under the Pfaffian 
then has two exact zero modes. The corresponding phasefactor then averages to 
zero exactly. 

5.2 D = 2 specialties 

Fermions in two Euclidean or one space dimension are simpler and in a way untyp- 
ical for the true problem related to the oscillating phase In the Euclidean field 
theory formulation this is seen by the phases from fermion loops and from spin 
'essentially canceling' in D = 2 (only). In our realization this is seen as follows. 
Minus signs appear only at two types of corners, namely (— 0| — 1) or ( — 1| — 0). 
By drawing closed loops with the intersection properties relevant here on a planar 
torus it is not difficult so see that 




(56) 



4 £) _ («^ i£) (k))) 



(57) 



z { p ((5 u , v ^'Kk))) 
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• loops winding around the torus in one or both dimensions receive an even 
number of such minus signs 

• loops that close trivially and do not wind around the torus receive an odd 
number of minus signs from spin phases. 

In [TT] a more detailed discussion of this and some illuminating figures with exam- 
ples can be found. Winding around the torus can be read off from the crossing of 
'boundary links' (1381) . Thus the result for each closed loop on the two-dimensional 
torus can be summarized in our notation as 

= L e ' (A,) = (0 ' 0) (B = 2only). (58) 

Negative signs — remember that the Fermi loop sign has been included in cf) - 
only come from topologically nontrivial loops that cannot be contracted to the 
trivial loop by series of plaquette moves [TUJ , [TT] . The total phase for u = v 
configurations can now be given as 

$( £ )(A;) = (_i)^+^,(o,o)+i ^ {oiu = l)) Q(k;u,u) = 1]. (59) 

It depends on k only via the topology variable e M . Using Fourier transformation 
on Z(2) E e (~l) e ^ e_e = 4<5 e ,e'] one may show the identity 

1 = z{e)1> (£) {k), z(e) = I(-l)*«.(o,o) (60) 

£ 

for this case. This in turn implies for the average monomer density (with no 
phases) 

the exact result (for free fermions) 
with the partition function 

Z = Y,<e)Z^ (63) 

e 

summed over boundary conditions with amplitudes z(e). Similarly from fl57|) we 
may deduce now 

tit = («^ £) (k))) m) 

Z ((5 U , V )) • 1 ] 
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6 Numerical experiments 



In this section we test our simulation method for the case of free fermions on 
various lattice sizes and (constant) m values in both D = 2 and 3. While of course 
no numerical simulations are really needed here, we nonetheless think that for 
our loop-gas representation this is not an untypical case also for later interacting 
applications. Thus the advance knowledge of the results here is just an advantage 
for precision testing. 

To extract fermionic quantities from simulations we must have the phase $^ (k) 
available for each sampled configuration for the desired boundary conditions. For 
each admissible configuration with Q(k; u, v) — 1 it can be constructed by tracing 
the string and all loops at the cost of order V operations. It is however easier to 
update its value together with the configurations. In fact this is even necessary 
to measure efficiently between microsteps as discussed in [TI]. To that end we 
assume & £ \k) to be known for the start configuration. We always took the trivial 
u — v = 0, ki = with & £ \k) = 1. Then, whenever an update proposal of type 
l v of the PS algorithm is accepted, we change 

$ (£) (A;) -> $ (e) (A/) = & £ \k) x ^- x r)((vv'),£) (65) 

\q\ 

and similarly for l u . Here q is given in tableland the additional factor 

/ / # \ \ f —1 if (vv 1 ) is au — boundary link and e u = 1 
V((w>),e) = { +1 elge " (66) 

takes into account the boundary conditions (see ( |38i) ). Needless to say, one may 
also keep track of $^ for several boundary conditions in the same run, as the 
updates do not depend on them. 

The spin matrix is easy to construct at any time from the single dimers adjacent 
to u and v ^ u and is trivial for the coinciding case. In practice we measure 
correlations contracted with some Dirac matrix T which leads to 

(Uil,M'Hk)tv[S(k;u,v)T])) 

-<«»)reW> = /K*)^ {{Svmm < 67 > 

We want to further specialize to zero spatial momentum as discussed in appendix 

m 

k T (z ) (mm*)) oc p(z ) ^_ V0!z y £ \k) tr^; u, v)T})) . (68) 

We took p to only depend on time and dropped the denominator. For symmetry 
reasons only T = 1, 7 (also labeled as S, V) leading to the scalar and vector correla- 
tions k$,ky are nontrivial. During the simulation we simply add the corresponding 



18 



amplitudes into bins for each separation zq and then end up with correlations mea- 
sured for all distances. Pre-computed tables are heavily used to speed up and they 
lead to a very simple code. In the simulations reported below we have observed 
Metropolis acceptance rates close to 50% for D = 2 and 30% for D = 3. This 
is close to the amplitude change by 1/(D + m) when the worm 'eats' a monomer 
which it has to do to grow. In equilibrium also the other processes are important, 
but this one seems to set the scale. 

All error estimates below are derived with the method and tools detailed in 
|20j . In particular the definition of integrated autocorrelation times r int employed 
here can be found, see also remarks in [13]. Due to time series of length 10 6 and 
more the convolution step in UWerr, eq. (31) in [20], became a bit slow for online 
data analysis. We therefore tailored a special version UWerr _f ft which accelerates 
this step by using the fast Fourier transform. It is available on the web under 
www . physik . hu-berlin . de/com/ ALPHAsof t. 




Figure 3: Two typical configurations on a 64 2 lattice at criticality (m = 0, left 
panel) and with correlation length 64/5 (m = 0.0812 . . . ., right panel). The string 
a is given by the fat (red line), the other lines are loops Xj. Readers are asked to 
identify left-right and top-bottom edges in their mind. 



6.1 D = 2, physically large volume 

For tests in this subsection we chose a mass such that ujL = 5 holds with the pole 
mass uj = ln(l + m). The zero momentum timeslice correlations (1681) then fall 
off exactly with exp(— uxq), modified to cosh or sinh due to time periodicity, see 
appendix |A] A typical configuration is visualized by the right picture in figure 
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El We see that for fermions, in contrast to the Ising model, there really is a 
unique 'worm', which moves by the updates^. In all simulations the complete zero 
momentum two-point functions at all separations were consistent within errors 
with the exact results. We routinely checked plots of the deviation in units of 
the estimated error against Xq which are order one, occasionally straddling ±2. 
In addition diagnostic quantities like (joTj) and (16"4"|) were monitored. Because in 
a large box (compared to the inverse mass) few configurations wind around the 
torus we find no significant difference between periodic and antiperiodic boundary 
conditions. In the example to follow we measured (|64p and obtained 

-!L- = 0.9744(5), -!L- = 0.9745(5) 
in agreement with the exact answer. 



1.01 



o 



2 > 1.005 



E 

i > 



0.995 



0.99 





Figure 4: Correlation function ky and the effective mass derived from it. Errorbars 
are one sigma high. 

As an impression for the reader we show in figure H] results for the vector 
correlation ky on a lattice L = 64, T = 2L with = (1, 0) after accumulating 10 7 

6 The poor unoriented Majorana worm has however no distinction between head and tail! 
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iterations (steps per site). The correlation length u 1 is hence about 13 lattice 
spacings. We have chosen the bias p as 



p(t) oc cosh (T/2 - t)} 



(69) 



which leads to a population of timeslices ((S Uo ~v ,t)) approximately flat in t. We 
refer the reader to the discussion in [14J, which can be taken over essentially 
without change. The upper panel shows the correlation ky(xo) itself normalized 
by its exact value. The growth of the errors from left to right is due to an increase 
of the integrated correlation time v int from about 0.6 iterations at short distances 
to about 15 iterations at xq = T/2. In the lower panel we give the effective mass 
as a function of distance by matching subsequent timeslices to 



Here errors appear (apart from xq very close to T/2) to be independent of the 
separation in agreement with the observed autocorrelations r int ~ 0.5 for all xq. 
The longer autocorrelations observed in ky apparently cancel in the ratio. From 
the fluctuations in figure H] we conclude qualitatively that statistical fluctuations 
at neighboring time separations are strongly correlated in ky, but much less so 
in m e fj. In a run with p = 1 the growth of r int for ky does not occur. Its error 
however grows in a similar way due to the larger variance coming from fewer data 
at large separation u — v (fewer 'long worms') when no bias p is applied. The more 
interesting effective mass is more accurate with the bias used for the figure. 

6.2 D = 2, physically small volume 

We now simulate at the critical point which in the free case is known to occur at 
m — 0. We remind that due to the infrared regulator given by the small inverse 
temperature T with antiperiodic boundary conditions e = 1 this is physically 
well-defined. Such applications are of interest in interacting theories to study 
nonperturbative renormalization using the universal finite volume continuum limit. 
To this end we report measurements of ks{T/A) and ky{T/2) with e = (1,0). 
Further motivation for the study of these objects can be found in appendix IA1 and 
refs. [TB], P2]. In table [3] we compile our results from performing 10 6 iterations at 
each of the lattice sizes. Again r int are given in iterations. 

Here the topology and the sign $ (1,0) fluctuate, but we can achieve a percent 
accuracy with the given statistics, which could be enlarged. 




< m cS = m cS (x + 1/2). 



(70) 



21 



T = L 


n b yj- i 


1 mi, Kg [1 /4J 


kv(T/2) 

ivy y-L i j 


7";^+ i. It 1 a\ 
' mt, Ky {1 /4J 


16 


0.010(7) 


1.85(8) 


0.997(10) 


1.11(2) 


32 


-0.011(15) 


3.6(4) 


1.007(16) 


1.43(4) 


64 


0.011(25) 


5.0(5) 


0.943(29) 


2.68(14) 



Table 3: Results at the critical point m — 0. The exact values for all L are 
ks{T/A) = due to chiral symmetry and &v(T/2) = 1 corresponding to canonical 
field normalization. 




Figure 5: A typical configuration on 24 x 12 2 with correlation length 2 (m = 
0.6487...). 

6.3 D = 3, sign problem 

It is trivial to adapt the code from two to three dimensions. On small lattices 
T, L = 4, 6 we performed similar validation tests as before with completely accurate 
and satisfactory results. Note that the formulae for k$, ky in appendix [A] are 
equally valid for D = 2, 3. It turns out, however, that now very suddenly as the 
volume is increased or the mass is lowered the sign fluctuations abruptly become 
so violent that no signal is left in (1671) and also for ratios of correlations as in (I70p 
all estimates yield '0/0' within errors: the sign problem. For large enough mass 
loops remain small and predominantly planar. Such loops are as in two dimensions 
with phase one. For a demonstration we show in table H] results for two cases just 
before trouble strikes. In figure [5] the last configuration of our run at uL = 6 
is shown. For uL = 4 and the same lattice size no meaningful results can be 
obtained anymore. Integrated autocorrelations times were close to 1/2 for all the 
quantities studied. Although the observables in the dimer ensemble are complex, 
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the averages of the imaginary parts vanish within errors as they have to, since they 
are parity odd. This was first checked and then used before forming quotients. A 
bias was not used here, p = 1. 





«<W» 


«*„,„*™» 
«<W» 


m cfr L(3.5) 


M/uj 


6 


0.1344(5) 


0.1343(5) 


6.04(5) 


0.51 


5 


0.0168(6) 


0.0164(6) 


4.2(10) 


0.27 



Table 4: Simulation results from 10 7 iterations on 24 x 12 2 lattices. 

In the last column we report a mass M which was extracted from the time 
slice correlation ((5 Uo - VOj t)) without phase factors and using the time periodic 5- 
function. It also shows a mass-plateau and, as we see, it is smaller than the physical 
mass. Due to interference effects the spinor correlation decays faster than this 
'geometric' one. From this observation one could think that the loops and strings 
have the 'wrong' size. We tried to generate them with a modified mass parameter 
m + Am such that M « uj is achieved and then reweighted the observables to 
the true mass, which is easy if the total monomer number is available. We find 
however that this simple idea does not improve the sign problem. 

7 Conclusions 

We have formulated the 'worm' algorithm of Prokof'ev and Svistunov for lattice 
fermions of the Wilson type with Wilson parameter r = 1. As in the Ising model 
it estimates stochastically by the Monte Carlo method the untruncated hopping 
parameter expansion of the partition function together with the graphs needed for 
the full two-point function, which can thus be computed. The PS algorithm very 
naturally lends itself to easily keep track of all phase factors and spin matrices that 
appear in the expansion. In two space time dimensions the contributions of all 
graphs are positive up to finite size effects and simulations are similarly efficient 
as in the Ising model. The all order hopping expansion is also worked out and 
numerically tested in three dimensions. Here the weights of fermion loops acquire 
complex phase factors and, for small mass and/or large volume lead to numerically 
uncontrollable fluctuations. The very sharp borderline was found around correla- 
tion length two for a 24 x 12 2 lattice. In particular, the continuum limit cannot 
be approached. Clearly here the method has to be complemented for instance by 
an improved estimator which sums some part of the contributions analytically as 
in cluster methods. No such method is known at present for the system at hand. 

For two dimensional fermions we now plan to add the interaction of the Gross- 
Neveu model. For the O(iV) invariant model the Majorana fermion discussed 
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here then has to be replicated iV-fold. For each multi-dimer configuration there 
are now between K{x) = and K{x) = N monomers at each site. The four 
fermion interaction can be enforced by integrating over the common external field 
m{x) site by site with the appropriate Gaussian weight yielding a i^(a;)-dependent 
total weight as already discussed in [UJ. In this way a coupling between the N 
'flavors' arises. The worm head and tail u and v now refer to one of the flavors. 
When u — v is reached with random re-location also a new flavor-index is chosen 
randomly. Further details still have to be worked out. 
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A Timeslice correlations of free fermions 

We prefer to re-introduce the lattice spacing a for this appendix. 

For our numerical test we need to know the free Wilson propagator (fT3|) in the 
time-momentum basis 

G(x , p k ) =a D - l J2 G(x)e-^ x » (71) 

where Latin vector indices k = 1,...,D — 1 refer to spatial components only. 
Spatial momenta are quantized as demanded by the spatial boundary conditions 
e k and sizes L k . From the Fourier expansion of G it is easy to see that G is given 
in terms of a function f(xo) by 

G(t,p k ) = {-iap klk + 1 + aM)f(t) - 1(1 + 7b)/(* + 1) - ^(1 - 7o)/(* - 1) (72) 

with ap^ = sin^apfj), M = m + ^p k p k , ap^ = 2 sm(ap ll /2). The function itself 
reads 

1 ^ e ipot 

f{t) = ^Tl^pi+p k p k+ (M+ipir (73) 

We find it useful to first evaluate / in the limit of infinite time extent T = Lq — > oo. 
We demand 1 + aM > and introduce 



sinh(au;/2) = — 



1 + aM 



au^O. (74) 
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With z = e iapo we now arrive at the contour integral 

/oo(t) = TT^M jtoriiz- e^)(z - e— ) (?5) 

and obtain 

f~>® = on 7 m • I-, e "^'- (76) 
2(1 + aiWj sinn au; 

The solution that is (anti)periodic over finite time now follows from 

oo 

f(t)= foo(t + nT)(-ir° (77) 

n=— oo 

and reads explicitly for ^ t < T 

1 / cosh[u;(T/2 - 01/ sinh(^T/2) for £ = 

M; 2(l + aM)sinhacc> \ sinh[^(T/2 - 1)}/ cosh(cuT/2) for e = 1 ' K ' 

For < t < T this yields 



rA/ M-aw 2 /2 f 

tr[G(g °' Pt)1= (l + a M)w X \ 



M - aw 2 /2 ^ / cosh[cu(T/2 - *)]/ sinh(a;T/2) for e = 
sinh[cj(T/2-t)]/cosh(cuT/2) for e = 1 

(79) 

with au> = 2 sinh(aa;/2), acu = sinh(aa;), and 



trhf _ n, _ 1 x / sinh[c(T/2 - 0]/ sinh(cT/2) for e = 

tr[7 G(a;o,Pfe)J - Y+~^M \ cosh[u;(T/2 - *)]/cosh(wT/2) for e = 1 " (80) 

At zero time separation tr[y o G(0,pk)] = vanishes and 

tr[G(0,p fc )] " T^m + (1 + aMp X \ tanh[cT/2] for e = 1 (81) 
Exact results for the zero momentum (for = 0) correlations 

k s (x ) = -a ' 1 £<£(0)£(a;)) = tr[G(x , 0)] (82) 

and 

M*o) = -a ' 1 £<£(0)7o£(aO> = tr[ 7o G(x , 0)] (83) 

follow from the above formulas. 

Incidentally these free field or tree level results form an example how Symanzik 

improvement for Wilson fermions works. If we multiply the left hand sides of 
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( 1751) and ( l80l) by the wavefunction improvement factor 1 + am and eliminate on 
the right hand side m by the pole mass m po i e = a -1 ln(l + am), then the finite 
size continuum limit (T/a — > oo at fixed t/T,m po \ e T,pkT) is reached at a rate 
proportional to (a/T) 2 . For p k = 0, where u> = m pole holds, lattice artifacts are 
even completely eliminated. 

In the correlations above one may also safely take the limit m — > at finite 
T for antiperiodic boundary conditions Eq — 1. This setup may be interpreted as 
a finite temperature which supplies an infrared regulator, similar as but simpler 
than in the Schrodinger functional [21] . In this limit ky becomes constant while 
ks vanishes. In D = 2 the latter may be viewed as a consequence of the discrete 
chiral symmetry £ — > «7o7i£ under which ks is odd. On the lattice with Wilson 
fermions this symmetry must emerge in the continuum limit at the critical mass, 
zero in the free case. Due to the (accidental) complete improvement it here holds 
already on the lattice at all nonzero separations. 

As discussed in [TO] ks, ky may be used to formulate massless nonperturbative 
renormalization conditions in the interacting Gross-Neveu model. The field can 
be normalized by keeping for instance ky(T/2) = 1 while the critical mass may 
be determined by tuning to ks(T/4) = 0. Note that ks{T/2) vanishes due to time 
reflection invariance combined with antiperiodicity for all mass values and is not 
suitable to fix the chiral point. 

B Spin phase: continuum and lattice 

We now consider an alternative way to compute the loop phase ( 1231) 

W (X) = - tr[P(mi)P(m 2 ) ■ ■ ■ P{m\\\)) (84) 

from the sequence of lattice unit vectors m,. An earlier discussion of this quantity 
is found in [22]. As a technical simplification we assume for a while that straight 
sections are absent, mj+i ■ mi — 0. We can rotate one projector into the next, for 
example 

P(mi)P(rri2) = R(mi,m2)P(m2)R(mi,m2}~ 1 P(m2) 

with the tc/2 spinor rotation matrix 

R(m 1 ,m 2 )=e* [vll ' vta] (86) 
in the usual notation yd = m M 7 M . By iterating this step we arrive at 

- w(X) = 2-(l Ah1 )/ 2 tT[R loop R(m u m w )P(m w )] (87) 



—=R(mi, m 2 )P(m 2 ) 
v2 



(85) 
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with 

-Rioop = R{m 1 ,m 2 )R(m 2 ,m 3 ) ■ ■ ■ R(m\x\-i, m\x\)R(m\\\, mi). (88) 

As the loop closes this total rotation must bring nil back to itself, it belongs to 
the stability group of this direction. 

The rotations introduced in ( 1861) can be characterized by mapping a cube onto 
itself by rotations by n/2 through any lattice plane. Products of such operations 
do not preserve the form (1861) but are general members of the cubic group. The 
number of elements is given by twice the number of orientations of a D-dimensional 
cube along the lattice axes (without reflections) that is 2 D D ! in D dimensions. The 
twofold multiplicity is because ±i? in the spinor representation corresponds to only 
one re-orientation of the cube as usual with spinor rotations. 

In D = 2 the cubic group is abelian and has 8 elements. Here -Ri 00 p adds up the 
half-angles of the rotations that a test-particle taken around the loop experiences. 
For a non-intersecting closed loop we get -Ri 00 p = — 1 for a 2ir rotation. Then we 
find in total 

w(\) = (_)-+i 2 -cW/ 2 ( 89 ) 

with the number of complete rotations v and the number of corners along the loop 
C(A). This form of the result now also covers the case where m^+i = rrii can occur, 
a straight section not counted in C(A). 

In D = 3 the cubic group is non-abelian and has 48 elements. Now -Ri 00 p can 
only be a rotation around the m\ axis, 

R loop = ±e*W (90) 

with a a multiple of n/2. Then one easily evaluates 

1 _. a 

ti[Ri oop R(m u m lxl )P(m lxl )] = -^=e (91) 

The sign ambiguity of the spinor transformation (2 to 1 covering) has been ab- 
sorbed into the a and the phase is in Z(8). For m\ = m\\\ again the factor l/y/2 
has to be dropped. 

It is interesting to note that a strictly analogous construction can be made 
in the continuum. Let us assume a closed curve j(t),t G [0,1]. For a regular 
parametrization we have 7^0 and form m(t) = 7/I7I, a mapping from Si to 
Sd-i- We now move m(t) along the loop by infinitesimal rotations in spinor form 

pi = pi pi 2 = -[pi pi, pi] (92) 

using m ■ m — 0. The solution of 

j t R(t) = ^pipiR(t), R(0) = l (93) 
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is formally given by the path-ordered integral 

R(t) = V exp j 1 jT drfA{r)fA{r) } . (94) 

The direction m is parallel transported according to 

p/i(t) = R(t)^)R-\t). (95) 

At t — 1 we must recover = pi(0) and thus R(l) must induce a rotation 

that has m(0) as a fixed point. In D = 2 there is no nontrivial such rotation, but 
since both R — 1 and i? = — 1 induce the trivial rotation of m(0), it is possible to 
find .R(l) = —1. Indeed the latter result emerges for a loop with an odd number 
of windings. In D = 3 R(l) can induce an arbitrary rotation around the m(0) 
axis, R(l) = exp[iapi(0)/2\. Thus there is a phase, now a general element of U(l), 
associated with a closed loop. It does not depend on where the loop is cut by the 
parametrization: going to another point, R(l) just gets similarity-transformed. 
It is a property just of the loop itself. It seems very likely that this is a known 
concept in differential geometry, however no reference is known to the author. 
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